library(tidyverse)
library(forcats)
library(dplyr)
library(foreign)
library(plotly)
library(survival)
library(flexsurv)  
library(aod)       
library(haven)
library(magrittr)
library(survminer)
library(lmtest)
library(sandwich)
library(stargazer)
library(MASS)
library(countrycode)
library(plm)
library(knitr)
library(pscl)
library(mediation)

## load dataset

westedu5 = read.csv("westedu_chaid_oct24.csv")

##main models

westedu5$cinc1 = westedu5$cinc*1000
westedu5$westedu2 = as.factor(westedu5$westedu2)

m1 =  lm(log.total.volume2 ~ westedu2
         +lag.gdp+ln_open+cinc1
         +polity2+Accountability
         +CNTS+Civil_War
         +forbfcol+gender
         +us_ally+unsc
         +entry
         +factor(leveledu)
         +businesspeople+teacher+lawyers+polactive
         +bluecollar+scientist+service+religious+milservice
         +ps.pa+economics+ossa+others
         +factor(ccode)+factor(year),
         data=westedu5)
m1.results = coeftest(m1, vcov. = vcovCL(m1, cluster = ~ccode+year, type = "HC0"))



m2 =  lm(log.oda.volume2 ~ westedu2
         +lag.gdp+ln_open+cinc1
         +polity2+Accountability
         +CNTS+Civil_War
         +forbfcol+gender
         +us_ally+unsc
         +entry
         +factor(leveledu)
         +businesspeople+teacher+lawyers+polactive
         +bluecollar+scientist+service+religious+milservice
         +ps.pa+economics+ossa+others
         
         +factor(ccode)+factor(year),
         data=westedu5)
m2.results = coeftest(m2, vcov. = vcovCL(m2, cluster = ~ccode+year, type = "HC0"))

m3 =  lm(log.oof.volume2 ~ westedu2
         +lag.gdp+ln_open+cinc1
         +polity2+Accountability
         +CNTS+Civil_War
         +forbfcol+gender
         +us_ally+unsc
         +entry
         +factor(leveledu)
         +businesspeople+teacher+lawyers+polactive
         +bluecollar+scientist+service+religious+milservice
         +ps.pa+economics+ossa+others
         
         +factor(ccode)+factor(year),
         data=westedu5)
m3.results = coeftest(m3, vcov. = vcovCL(m3, cluster = ~ccode+year, type = "HC0"))


stargazer(m1.results, m2.results, m3.results,
          covariate.labels = c("Western education", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
                              ))
stargazer(m1,m2,m3)

## negative binomial models

m4 =  glm.nb(total.count ~ westedu2
             +lag.gdp+ln_open+cinc1
             +polity2+Accountability
             +CNTS+Civil_War
             +forbfcol+gender
             +us_ally+unsc
             +entry
             +factor(leveledu)
             +businesspeople+teacher+lawyers+polactive
             +bluecollar+scientist+service+religious+milservice
             +ps.pa+economics+ossa+others
             +factor(ccode)+factor(year),
         data=westedu5)
m4.results = coeftest(m4, vcov. = vcovCL(m4, cluster = ~ccode+year, type = "HC0"))

m5 =  glm.nb(oda.count ~ westedu2
             +lag.gdp+ln_open+cinc1
             +polity2+Accountability
             +CNTS+Civil_War
             +forbfcol+gender
             +us_ally+unsc
             +entry
             +factor(leveledu)
             +businesspeople+teacher+lawyers+polactive
             +bluecollar+scientist+service+religious+milservice
             +ps.pa+economics+ossa+others
             +factor(ccode)+factor(year),
         data=westedu5)
m5.results = coeftest(m5, vcov. = vcovCL(m5, cluster = ~ccode+year, type = "HC0"))

m6 =  glm.nb(oof.count ~ westedu2
             +lag.gdp+ln_open+cinc1
             +polity2+Accountability
             +CNTS+Civil_War
             +forbfcol+gender
             +us_ally+unsc
             +entry
             +factor(leveledu)
             +businesspeople+teacher+lawyers+polactive
             +bluecollar+scientist+service+religious+milservice
             +ps.pa+economics+ossa+others
             +factor(ccode)+factor(year),
         data=westedu5)
m6.results = coeftest(m6, vcov. = vcovCL(m6, cluster = ~ccode+year, type = "HC0"))

stargazer(m4.results, m5.results, m6.results,
          covariate.labels = c("Western education", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
          ))
stargazer(m4,m5,m6)

## Count calculation: oda
model.data <- augment(m5) %>% 
  mutate(index = 1:n()) 
model.data = rename(model.data, year='factor(year)')
model.data = rename(model.data, ccode='factor(ccode)')
model.data = rename(model.data, leveledu='factor(leveledu)')

m5 = glm.nb(oda.count ~ westedu2
            +lag.gdp+ln_open+cinc1
            +polity2+Accountability
            +CNTS+Civil_War
            +forbfcol+gender
            +us_ally+unsc
            +entry
            +factor(leveledu)
            +businesspeople+teacher+lawyers+polactive
            +bluecollar+scientist+service+religious+milservice
            +ps.pa+economics+ossa+others
            +ccode+year,
            data=model.data)



education <- c('0','1')

qi <- sapply(education, FUN=function(x){
  mean(predict(m5, type="response", newdata=mutate(model.data, westedu2=x)), na.rm=TRUE)
})
qi <- data.frame(westedu = education, count = qi)

## Count calculation: oof
model.data <- augment(m6) %>% 
  mutate(index = 1:n()) 
model.data = rename(model.data, year='factor(year)')
model.data = rename(model.data, ccode='factor(ccode)')
model.data = rename(model.data, leveledu='factor(leveledu)')

m8 = glm.nb(oof.count ~ westedu2
            +lag.gdp+ln_open+cinc1
            +polity2+Accountability
            +CNTS+Civil_War
            +forbfcol+gender
            +us_ally+unsc
            +entry
            +factor(leveledu)
            +businesspeople+teacher+lawyers+polactive
            +bluecollar+scientist+service+religious+milservice
            +ps.pa+economics+ossa+others
            +ccode+year,
            data=model.data)



education <- c('0','1')

qi <- sapply(education, FUN=function(x){
  mean(predict(m6, type="response", newdata=mutate(model.data, westedu2=x)), na.rm=TRUE)
})
qi <- data.frame(westedu = education, count = qi)

## no control

m7 =  lm(log.total.volume2 ~ westedu2,
         data=westedu5)
summary(m7)

m8 =  lm(log.oda.volume2 ~ westedu2,
         data=westedu5)
summary(m8)

m9 =  lm(log.oof.volume2 ~ westedu2,
         data=westedu5)
summary(m9)

m10 =  glm.nb(total.count ~ westedu2,
         data=westedu5)
summary(m10)

m11 =  glm.nb(oda.count ~ westedu2,
         data=westedu5)
summary(m11)

m12 =  glm.nb(oof.count ~ westedu2,
         data=westedu5)
summary(m12)

stargazer(m7, m8, m9, m10, m11, m12)

## Western aid

westedu51 =  read_csv("westedu_westaid_oct24.csv")
westedu51$log.west.aid = log(westedu51$west.aid+1)
westedu51$log.west.aid1 = westedu51$log.west.aid
westedu51$log.west.aid1[is.na(westedu51$log.west.aid1)] = 0

westedu51$cinc1 = westedu51$cinc*1000

m13 =  lm(log.west.aid ~ westedu2
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
         +factor(ccode)+factor(year),
         data=westedu51)
m13.results = coeftest(m13, vcov. = vcovCL(m13, cluster = ~ccode+year, type = "HC0"))

stargazer(m13.results,
          covariate.labels = c("Western education", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
          ))
stargazer(m13)

## total aid

m14 =  lm(log.aid.amount ~ westedu2
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
          data=westedu51)
m14.results = coeftest(m14, vcov. = vcovCL(m14, cluster = ~ccode+year, type = "HC0"))

stargazer(m14.results,
          covariate.labels = c("Western education", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
          ))
stargazer(m14)

## Accountability

m15 =  lm(Accountability ~ westedu2
         +ln_open+lag.gdp
         +factor(leveledu)
         +CNTS+forbfcol
         +Civil_War+gender
         +military+businesspeople+teacher+lawyers+polactive
         +bluecollar+scientist
         +polity2
         +ps.pa+economics+ossa+others
         +entry+cinc+us_ally+unsc
         +factor(ccode)+factor(year),
         data=westedu5)
summary(m15)
coeftest(m15, vcov. = vcovCL(m15, cluster = ~ccode+year, type = "HC0"))

m16 =  lm(log.total.volume2 ~ westedu2+Accountability
         +ln_open+lag.gdp
         +factor(leveledu)
         +CNTS+forbfcol
         +Civil_War+gender
         +military+businesspeople+teacher+lawyers+polactive
         +bluecollar+scientist
         +polity2
         +ps.pa+economics+ossa+others
         +entry+cinc+us_ally+unsc
         +factor(ccode)+factor(year),
         data=westedu5)
summary(m16)

stargazer(m15, m16)
mediation2 <- mediate(m15, m16, sims=1000, 
                      treat="westedu2", mediator="Accountability")
summary(mediation2)

## type of education

westedu5$edutype = 0
westedu5$edutype = ifelse(westedu5$degree=='College', 2, westedu5$edutype)
westedu5$edutype = ifelse(westedu5$degree=='College and postgraduate', 2, westedu5$edutype)
westedu5$edutype = ifelse(westedu5$degree=='High school and college', 2, westedu5$edutype)
westedu5$edutype = ifelse(westedu5$degree=='Master', 1, westedu5$edutype)
westedu5$edutype = ifelse(westedu5$degree=='Master and PhD', 1, westedu5$edutype)
westedu5$edutype = ifelse(westedu5$degree=='Phd', 1, westedu5$edutype)

westedu5$edutype = as.factor(westedu5$edutype)



m17 =  lm(log.total.volume2 ~ edutype
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
         data=westedu5)
m17.results = coeftest(m17, vcov. = vcovCL(m17, cluster = ~ccode+year, type = "HC0"))

stargazer(m17.results,
          covariate.labels = c("Western education: only Master/PdD", "Western education: College", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
          ))
stargazer(m17)

## Africa subsample

westedu4 = subset(westedu5, region1==11)

m18 =  lm(log.total.volume2 ~ westedu2
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
          data=westedu4)
m18.results = coeftest(m18, vcov. = vcovCL(m18, cluster = ~ccode+year, type = "HC0"))

stargazer(m18.results,
          covariate.labels = c("Western education", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
          ))
stargazer(m18)

## transparency

trans = read_csv("transparency.csv")


TransparencyData$ccode = countrycode(TransparencyData$Country, origin = "country.name", destination = "cown")
disaster = mutate(disaster, ccode=replace(ccode, Country=="Serbia", 345)) 
disaster = mutate(disaster, ccode=replace(ccode, Country=="Yemen", 679)) 

trans = rename(trans, year=Year, trans=`HRV Transparency index`)
TransparencyData = dplyr::select(TransparencyData, ccode, year, transparencyindex)

westedu5 = left_join(westedu5, TransparencyData)

m19 =  lm(transparencyindex ~ westedu2
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
          data=westedu5)
coeftest(m19, vcov. = vcovCL(m19, cluster = ~ccode+year, type = "HC0"))
summary(m19)


m20 =  lm(log.total.volume2 ~ westedu2+transparencyindex
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
          data=westedu5)
summary(m20)

stargazer(m19,m20)

mediation3 <- mediate(m19, m20, sims=1000, 
                      treat="westedu2", mediator="transparencyindex")
summary(mediation3)


## interaction

m21 =  lm(log.total.volume2 ~ westedu2*businesspeople
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
          data=westedu5)

m22 =  lm(log.total.volume2 ~ westedu2*lawyers
          +lag.gdp+ln_open+cinc1
          +polity2+Accountability
          +CNTS+Civil_War
          +forbfcol+gender
          +us_ally+unsc
          +entry
          +factor(leveledu)
          +businesspeople+teacher+lawyers+polactive
          +bluecollar+scientist+service+religious+milservice
          +ps.pa+economics+ossa+others
          +factor(ccode)+factor(year),
          data=westedu5)

m21.results = coeftest(m21, vcov. = vcovCL(m21, cluster = ~ccode+year, type = "HC0"))
m22.results = coeftest(m22, vcov. = vcovCL(m22, cluster = ~ccode+year, type = "HC0"))

stargazer(m21.results, m22.results,
          covariate.labels = c("Western education", "Business", "GDP per capita", "Trade openness",
                               "National Capability", "Polity", "Accountability", "CNTS",
                               "Civil War", "Before college", "Gender", "US ally", "UN Security Council"
          ))
stargazer(m21, m22)

library(margins)
quantile(spend5$v2x_polyarchy, probs = c(0.1, 0.25, 0.5, 0.75, 0.9), na.rm=TRUE)
margins1 <- summary(margins(m20, variables = "westedu2", at = list(military = c(0,1))))

LOWER = margins1$AME + qnorm(.95) * c(-margins1$SE)
UPPER = margins1$AME + qnorm(.95) * c(margins1$SE)
margins = cbind(margins1, LOWER, UPPER)
margins["quantile"] = c(0.1, 0.25, 0.5, 0.75, 0.9)
g4 <- ggplot(margins, aes(x=factor(quantile), y=AME, group=1)) +
  geom_point() +
  geom_errorbar(aes(ymin=LOWER, ymax=UPPER), width=.2,
                position=position_dodge(0.05)) +
  geom_hline(yintercept = 0, lty=2) +
  xlab("Electoral democracy") +
  ylab("Marginal effect of domestic attack") 
g4

## westedu aid

westedu6 = aggregate(westedu5$total.volume1, by=list(westedu5$year, westedu5$westedu2), FUN=sum)

westedu6 = rename(westedu6, year=Group.1, Western_education=Group.2, total.amount=x)

westedu6$Western_education = as.factor(westedu6$Western_education)

g <- ggplot(westedu6, aes(x=year, y = total.amount, group=Western_education, color=Western_education)) +
  geom_line() +
  geom_point() +
  xlab("Year") +
  ylab("Total amount of Chinese OF")
g
